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We present a new procedure which enables to extract a scattering matrix S {I) as a 
complex function of angular momentum directly from the scattering data, without any a 
priori model assumptions implied. The key ingredient of the procedure is the evolutionary 
algorithm with diffused mutation which evolves the population of the scattering matrices, 
via their smooth deformations, from the primary arbitrary analytical S (l) shapes to the final 
ones giving high quality fits to the data. Due to the automatic monitoring of the scattering 
matrix derivatives, the final S (l) shapes are monotonic and do not have any distortions. 
For the ^®0-^^0 elastic scattering data at 350 MeV, we show the independence of the final 
results of the primary S (l) shapes. Contrary to the other approaches, our procedure provides 
an excellent fit by the S (l) shapes which support the "rainbow" interpretation of the data 
under analysis. 

PACS number(s): 24.10.Ht, 25.70.-z, 25. 70. Be 



I. INTRODUCTION 



S'-operator is a fundamental quantity of the scattering theory, which incorporates, by a general assumption, 
all possible information on any possible scattering process (including particle creation/destruction). In 
the case of an elastic scattering, the diagonal matrix elements of S'-operator in the angular momentum 
representation can be given in general form as 

S(0-?7(/)e2'^W, (1) 

where the S'-matrix modulus rj (l) and the scattering phase ip (l) are real smooth functions of I. The unitarity 
of the S'-matrix for the composite-particle-nucleus scattering in the presence of nuclear absorption requires 
that 7/ (Z) < 1, so we put 

77(0 = e-2^"«, (2) 

where the nuclear absorption phase 6a (l) must be a real smooth positive function of I. 

Since the colliding nuclei have electric charges, then the scattering phase tp (I) can be divided into two 
parts 

^{l)^Sr{l)+ac{l), (3) 

where the nuclear refraction phase 5r (l) and the Coulomb scattering phase ac (0 must be real smooth 
functions of I. 

From a general physics viewpoint, the only restrictions we may impose on the nuclear phases 6a.r (0 to 
be determined are their finite values at small I, total vanishing at sufficiently large I and smooth behavior 
in the intermediate region. The most natural and simple approximation for da (0 (or 77 (/)) and dr (l) is a 
monotonically descending (for 77 (l), ascending to unity) function which can be easily modelled with help of, 
say, the Fermi-step or Gauss functions. For the case of elastic heavy-ion scattering at intermediate energies 
(>20 MeV/nucleon), the S-matrix approaches of such a kind (see, e.g., [1], [2], [3], [4]) and the optical 
potential models which yield S (l) with such a behavior (see, e.g., [5], [6]) have appeared quite successful and 
argued for the so-called "rainbow" interpretation of the data. However, these models have not allowed the 
adequate description of all the features of the data measured. 

At the same time, in many cases the quality of fit can be improved when the phases Sr {l)aTe modified by 
the additional surface terms of different forms (see, e.g., [7], [8]). Such modifications, in general, make the 
S (Z)dependence nonmonotonic. Note also the S-matrix model with the additional derivative-like interior 
term in the absorption phase Sa (0 [9]- The nonmonotonic behavior of the described type is also inherent 
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in the scattering matrices found with help of the optical potentials which have both the standard Saxon- 
Woods forms and the ones with the additional surface terms (see, e.g., [6], [10], [11]). In spite of the 
nonmonotonicity of S (l) in these approaches, the mentioned above "rainbow" interpretation of the data is, 

nevertheless, preserved. 

Further substantial improvement in the quality of fit is achieved with help of the more flexible S (/) forms 
which allow the phases to behave nonmonotonically for all relevant I. Such a nonmonotonic behavior is 
provided by extending the standard (monotonic) 5-matrices with the scries of the pole- like terms (sec, e.g., 
[12]) or the proper (say, spline) basis functions (see, e.g., [13], [14], [15]). The similar behavior is inherent in 
the ^-matrices calculated from the optical potentials which have the additional derivative-like interior terms 
or have the more complicated forms obtained with use of the spline functions or the Fouricr-Besscl scries 
(see, e.g., [5], [16], [17], [18]). In spite of the excellence of the quality of fit provided in such approaches, the 
"rainbow" interpretation of the data appears to be no longer valid, which rises the problem of finding the 
physical meaning of the results obtained this way. 

Clearly, all the approaches mentioned above are more or less model-dependent because the functions used 
to model the phases (5„.,. (/) and the real and imaginary parts of optical potential V (r) and W (r) are more 
or less the properly parametrized analytical ones. Thus, the search spaces of all possible shapes for the 
5-matrix and the optical potential are strongly reduced and, consequently, the data analyses performed on 
such spaces can lead to the incorrect physical interpretation of the data. 

That is the reason why it would be highly desirable to have the procedure which could be able to extract 
the scattering matrix and/or the optical potential directly from the experimental data, without introduction 
of any bias towards some a priori "physically reasonable" model assumptions. The very first question this 
procedure must answer to is whether the nonmonotonic (or pole-like) structures and any other distortions, 
which appear in the 5'-matrix shapes obtained in the most successful approaches, are really necessary to 
reproduce the experimental data studied. This will help us to shed more light on the applicability of the 
"rainbow" interpretation to the heavy-ion collisions in the wide range of energies and mass numbers. 



II. MODEL-FREE DETERMINATION OF THE SCATTERING MATRIX 



To develope the desired procedure which determines S (1) directly from the data, data — > we 
need to solve the problem in its most explicit form where each value of Sa^r (0 is treated, generally, as an 
independent fitting parameter. This makes the problem parameter space highly dimensional and the choice 
of an appropriate search method crucial. Evolutionary (or genetic) algorithms (EA's) have been many times 
proven very efficient in dealing with very difficult physical problems (see, e.g., [19], [20], [21], [22]), so we 
have chosen E A as a key element of our procedure. Note that our algorithm resembles the so-called smooth 
genetic algorithm proposed in [23] . 

According to the general ideology of the EA implementation, we deal with the population of N indi- 
viduals. Each individual is the S'-matrix presented as the pair of the real-valued /max-dimensional vectors 
{Sa (0 , Sr (/)), I = 0,1, Zmax — 1- The fitness of each individual reflects the quality of data fitting provided 
by the individual's ^-matrix. By using the mutation operation the algorithm evolves the initial population 
of the badly fitted individuals to the final population of the highly fitted ones. 

Every iteration of our procedure contains the following steps. 

1. Generate the initial population of A'^ individuals. For each individual, the vectors 5a.r (0 arc filled with 
help of any monotonically descending function of I, the first derivative of which has the only one minimum. 
To be definite and to test the robustness of the procedure against various starting conditions, we choose the 
following five primary models for S (l). 

a). The 6-parameter model composed of two Fermi functions: 



2Si{l) = Qi f {l,li,di) , f{l,li,di) 



1 + exp 



I - I, 



a, r. 



b). The 4-parameter model composed of two Gauss functions 

di 



25i (0 = gi exp ( ) • 



c) . The 5-parameter Mclntyre model [1]: 

r; (0 = / (-/, -la, da) , 26r {I) = Or f [l, Ir, dr) . 

d) . The 6-parameter phenomenological model [4]: 

25, (0 = .9, [{21 + 1) d,F {I, h, d,) + d^i^^ ^.-)] 1/2 ^ 

F {I, k, di) = -/-I (l, k, di) In [1 - / (/, k, di)] , 
r)„ = 1 . r)„ = 2. 



(4) 



(5) 



(6) 



(7\ 



e). The 6-parameter model composed of two power- type functions: 



2-5* (0 = y^nfTg-' > 2- (8) 

The parameters gi, li, di, and arc positive. They are chosen for each individual and each model 
function at random within some intervals which are wide enough to produce substantially different shapes of 
the phases. Normally, all the individuals in a given population are initialized with one and the same function 
from the set a) - e). 

All the mentioned primary models for S (l) are "physically justified" , except for the case e) which has 
no physical background. Nevertheless, wc have included this purely mathematical case to see whether the 
procedure is able to find "physically meaningful" results under such a tough conditions. 

2. Evaluate the fitness of each individual in the population. The fitness function in our approach consists 
of two parts. The first one is associated with the quality of the shapes of 6a,r (0 while the second one 
accounts for the quality of the fitting of the experimental data. 

The requirements which the shapes of 6a^r (0 must meet in our approach are as follows. 

i) . The functions Sa,r (l) must be descending. 

ii) . The first derivatives of 6a,r (0 ™ust have only one minimum and no maxima. 

iii) . The second derivatives of Sa,r (0 must not have more than one minimum and one maximum. 

iv) . The third derivative of 5r (l) must not have more than one minimum and one maximum. 

v) . The logarithmic derivative of 6r (l) must be descending. 

The requirements i) - iii) ensure the absence of any distortions of the phase shapes, at least, up to 
the derivatives of the second order. The condition iv) is added because we want the defiection function 
9 (l) = 2d(p (l) /dl to have no shape distortions up to the same order of its derivatives. The condition v) 
provides for the permanent decrease of 5r {I) with the increase of I. The requirements i) - iv) are crucial for 
the shapes of 5a,r (0- Thus the penalties imposed upon the individual in the case of the violation of these 
requirements are fatal. The condition v) is not so strong and introduces only the ultimate bias towards the 
desired tail of (5,. [l). 

The quality of the fit of the calculated differential cross section to the experimentally measured one is 
assessed via the standard magnitude per data point. The calculations are made by using the expansion 
of the scattering amplitude into a series of Legendrc polynomials.. The elastic scattering differential cross 
section is equal to the squared modulus of this amplitude. 

It is often claimed that the amount of the large scattering angle data is insufficient to determine the 
scattering matrix and/or the optical potential in a unique way. Thus, we add several additional pseudo data 
points after the last actual ones, which follow the tendency of the cross section behavior (cf., e.g., [14]). Of 
course, this prescription can not be universal and must be used with care in the context of the data under 
study. The incorporation of the invented data points to the criterium can appear misleading for the 
fitting procedure, therefore, we use the penalty-free corridor around those points and apply the prescription 
only after the fitting to the actual data set has been accomplished. 

3. Let each individual in the population produce M offsprings. The replication is performed according to 
the transformation: 



log {5[ {I)) = log (<5, (0) + A, (0, 1) D {I, l„ 



i), 



(9) 



where 5i (l) and S'^ (l) are the parent's and offspring's S'-matrix phases, respectively, > is the mutation 
amplitude, Ni (0, 1) denotes a normally distributed one-dimensional random number with mean zero and one 
standard deviation, Im.i stands for the mutation point chosen randomly in the interval < ImA < ^max — 1, 
and dm.i > is the value characterizing the diffuseness of the mutation point. The diffusing function 
D{1, lm,i, dm,i) must be of the bell-like shape with the only maximum at I = lm,i and the fall-off tail 
around this point. To be definite and to ensure the proper localization of the consequences of the mutation 
we choose the diffusing function in the form: 



^ ij'i ^m,ij dm,i) — exp 



(10) 



The mutation amplitude Ai and the mutation diffuseness d^.i are the quantities automatically tuned 
within some intervals. The limits of these intervals, having the extremely large values at the beginning of 
the procedure, are smoothly decreased in the course of the run, and acquire the small values at the end. 
Such a schedule provides for both the removal of the features of the primary parameterizations a) - e) from 
the individual's S (l) and the fine tuning of the details of S (1). 

4. Evaluate the fitness values of all offsprings. Sort the offsprings in a descending order due to their 
fitnesses. Select best offsprings to form the new population. 

5. Go to step 3 or stop if the best fitness in the population is sufficiently high (the value is small 
enough). 



Evolutionary algorithms arc, generally, the global optimization technique which, however, can not guaran- 
tee that the optimum found is the global one. Therefore, it is necessary to run the procedure several times. 
Besides, there is no way to know in advance what will be the minimum value of the magnitude. Thus, 
it is instructive to monitor the dynamics of the best, worst and mean fitness values and the rms deviation 
from the mean fitness in the population during those several runs of the procedure. Such monitoring usually 
helps to localize the region of the potentially lowest values. 

III. SCATTERING MATRIX FOR THE ^^O-^^O ELASTIC SCATTERING AT 350 MEV 

We have applied our technique to analyze the well known test case of the ^^O-^^O elastic scattering at 
350 MeV, for which the approaches giving very good quality of fit predict the existence of the nonmonotonic 
structures in the 5-matrix (sec. e.g.. [12], [13]). 

In our calculations, bearing in mind that the collision energy is sufficiently high, we let ac (l) in Eq. (3) 
to be the quasiclassical phase of the point-charge scattering by the uniformly charged sphere (see, e.g., [3]) 
having the radius Rc = 0.95 x 2 x 16^/'^ [24]. The calculated differential cross sections have been symmetrized 
for the scattering of identical nuclei. Besides, the experimental errors are assumed to be equally weighted 
(10% error bars). 

Figures 1-5 show the results of our calculations with the primary models a) - e) for S (l), respectively. 
The x^ values for our fits to the data are 2.4 — 2.5. For each initial case, the results of five different runs 
of the procedure are presented to display the error bands within each of the primary S (/) models. Figure 6 
compiles five best results from Figs. 1-5 to illuminate their sensitivity to the details of the particular primary 
S (Z) model. Figure 7 demonstrates the consequences of the consideration of the invented data points in the 
region of large scattering angles. 

IV. DISCUSSION 

The evolutionary procedure of determining the scattering matrix, presented in this publication, is aimed 
to search for the globally optimal solution. But, being aware of the complexity of the problem under study 
and the fact that the actual number of fitting parameters (twice the number of angular momenta which is 
Zmax =120 in our test case) is substantially greater the actual number of data points (which is equal to 105 
in our test case), we do not expect to achieve it. Therefore, we consider the obtained results (Figs. 1-7) as 
very promising. 

First of all, we see that within every model used for the primary S (/) dependence, regardless the variety of 
their shapes, the moduli rj (1) and the nuclear refraction phases Sr (/), as well as the total deflection functions 
9 (l), obtained in different runs of the procedure, go close to each other (Figs. 1-5). The differences between 
them can sometimes be seen only in the enlarged or even logarithmic scale. The same observation can be 
made if one analyzes the compilation of the best results (Fig. 6), which points on their independence upon 
the initial conditions. At the same time, the nuclear phases 6r (0 deviate from each other in the region of 
large angular momenta. There the scattering matrix moduli r] {I) are very close to unity, which makes the 
contributions of the partial waves with these values of / to the scattering amplitude vanishingly small. To 
introduce the corresponding bias into the searching procedure, we probably need more precise experimental 
information in the region of small scattering angles. Nevertheless, we are able to conclude that, under the 
requirements i) - v) imposed upon the phases 6a,r (0> have managed to localize the region of the scattering 
matrix shapes which give the lowest values to the x^ magnitude. 

It is important to emphasize the remarkable fact that the application of the power-type function e) [Eq. 
(8)] as the primary S (l) parameterization, which has no proper physical meaning, does not produce any 
difficulties for our procedure to find the physically meaningful scattering matrix (Fig. 5) . From the formally 
mathematical viewpoint, the iterative application of the replication transformation (9) to the phases Sa^r (0 
is equivalent to the addition to the primary phases of the ultimately "infinite" sum of the diffusing functions 
(10) with various parameters and weights. Due to the special schedule of choosing and tuning the latter, 
the phase shapes are transformed almost adiabatically across the run of the procedure. As the result, the 
phases evolve to the equilibrium, with respect to the fitness, shapes which are free from any recollections 
about the particular models for the primary S {I) and the diffusing function. That is why, we believe, our 
procedure is actually a model-free one. 

Somewhat surprising seems the observation that the incorporation of the additional pseudo data points 
in the fitness function, which really forces the cross section to behave as desired, produces no noticeable 
corrections to the scattering matrix (Fig. 7) in the whole range of I. This is against the conventional way of 
thinking but can be just the feature of that particular data set under study. 

From the physics viewpoint, our results support the "rainbow" interpretation of the given data: the 
maximum in the differential cross section observed near 50° is identified as the primary nuclear rainbow. 
The nuclear rainbow angle which corresponds to the minimum of the defiection function 6 {I) acquires the 
values 9^ = 61 — 64°. At this noint, one might ask whether it is Dossible to imnrove the oualitv of fit, bearine 



in mind the number of fitting parameters. In fact, the answer is positive. If, from the very beginning, we 
abandon all the requirements i) - v) imposed on the shapes of 5a,r (0) ^^^^ the procedure becomes able to 
find the results with ~ 0.5 — 0.6. But the ^-matrices for these cases are nonmonotonic and substantially 
different from run to run, belonging to difFercnit local optima. The other way to search for the better quality 
of fit could be found in testing the stability of the monotonic shapes of the ^-matrices obtained in our 
study against the nonmonotonic transformations (9). Then it seems more probable to find the results which 
belong to the same local optimum or the nearest ones. Following this way, if and only if the substantial 
improvement in the quality of fit is accompanied by the repeated observations of the same equilibrium 
nonmonotonic structures in S{1), then the appearance of these structures should be admitted as necessary 
and the search for their physical interpretation urgent. 

The proposed recipe can also be useful in analyzing some important cases where the presence of nonmono- 
tonic (or even nonsmooth) structures in the ^-matrix seems to be justified. Namely, these are the cases 
where, for instance, the behavior of S {I) is resonance-dominated (see, e.g., [6], [27]), or it is important to 
account for the dynamic effects (parity dependence of the interactions between nuclei, elastic transfer, see, 
e.g., [28], [29]), or the interference effects condition the nonmonotonic scattering matrices (sec, e.g., [10]). In 
order to confirm the existence of the discussed equilibrium structures in the .S'-matrix, we need to extract 
them directly from the respective experimental data, using, for instance, our approach with the requirements 
i) - v) switched off cither from the very beginning or after the monotonic shape of the S'-matrix is obtained. 
If we succeed, then we must admit that the requirements i) - v) cannot apply to all cases. 

The evolutionary procedure presented in this publication has been initially devised to determine the 
scattering matrix in the angular momentum representation. Obviously, the similar approach can be used 
to develop the evolutionary procedure for the determination of radial dependence of a complex optical 
potential. With help of this procedure, the optical potential can be extracted directly from the experimental 
data: data — > V (r). Moreover, using the similar procedure, the scattering matrix produced by the optical 
potential can be fitted to the scattering matrix extracted directly from the data: data — > S (l) — > V (r). 
This means that the optical potential found this way will correspond to the scattering matrix extracted 
immediately from the data. Having unified these three search procedures data — > <S'(/), data — > ^ir), 
and data — > S (Z) — > V (r) into the one, we obtain a powerful tool for the deep theoretical investigation of 
the heavy-ion collisions at intermediate energies. 
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Figure Captions. 



FIG. 1. Five scattering matrices for the ^^O-^^O elastic scattering at 350 MeV, calculated by our procedure 
with the primary model a) for S (l) [Eq. (4)]. (a) Scattering matrix moduli rj (/). The inset shows the region 
of small momenta in the logarithmic scale, (b) Nuclear phases 6r (l). The inset shows the region of small 
momenta in the enlarged scale, (c) The same as (b) but in the logarithmic scale, (d) Deflection functions 
9 (l). The inset shows the vicinity of 8 (Z) minima in the enlarged scale. Solid curves correspond to the best 
quality of fit to the data = 2.4. 

FIG. 2. The same as FIG. 1 but with the primary model b) for S (l) [Eq. (5)]. 

FIG. 3. The same as FIG. 1 but with the primary model c) for S (l) [Eq. (6)]. 

FIG. 4. The same as FIG. 1 but with the primary model d) for S (/) [Eq. (7)]. 

FIG. 5. The same as FIG. 1 but with the primary model e) for S (l) [Eq. (8)]. 

FIG. 6. Five best results from Figs. 1-5. Notations are the same as in FIG. 1. 

FIG. 7. Two scattering matrices and differential cross sections for the ^^O-^^O elastic scattering at 350 

MeV, calculated by our procedure with the primary model b) for S (l) [Eq. (5)]. Solid (dashed) curves 
are the results of calculations with the invented data points taken (not taken) into account in the region 
of large scattering angles, (a) Scattering matrix moduli r] (l) and (b) nuclear phases Sr (1) in the region of 
small momenta, (c) 5,, (l) in the region of large momenta, (d) Deflection functions G (l) in the vicinity of the 
minima, (e) The differential cross sections (ratio to Rutherford). Experimental data are taken from Refs. 
[25] and [26]. Solid curves presenting S (1) correspond to the same ones shown on FIG. 2. 

Note to the Editor. Please, place each figure on a single journal page. 
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